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This  paper  develops  particle-resolved  simulations  to  predict  conductivity  within  porous  composite 
electrodes.  Hundreds  of  spherical  particles  (order  fractions  of  a  micron)  are  packed  randomly  into  a 
cubical  region  (order  of  a  few  microns),  using  two  alternative  packing  algorithms.  The  composite  struc¬ 
tures  include  both  ion-conducting  and  electron-conducting  particles.  The  particle  network  is  discretized 
using  tetrahedral  meshes  that  fully  resolve  the  interiors  of  the  particles  and  their  intersections.  Charge- 
conservation  equations  are  solved  to  predict  current  through  the  network.  These  simulations  are  used 
to  derive  the  effective  conductivities  that  are  required  for  macroscale  simulations  at  length  scales  much 
larger  than  the  particle  scale.  Because  the  microstructures  are  synthesized  via  random  particle  packing, 
multiple  realizations  are  needed  to  deliver  statistically  invariant  results.  The  results  show  that  a  few 
hundred  particles  with  a  few  hundred  realizations  is  sufficient.  Predicted  coordination  numbers,  perco¬ 
lation  probabilities,  and  three-phase-boundary  lengths  are  consistent  with  percolation  theory,  but  the 
predicted  effective  conductivities  are  significantly  smaller  than  those  predicted  with  conventional  per¬ 
colation  theory.  By  adjusting  the  Bruggeman  factor  from  the  conventional  value  of  1.5-3. 5  brings  the 
percolation-theory  prediction  for  effective  conductivity  in  line  with  the  fully  resolved  results. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  objective  of  this  paper  is  to  develop  microstructural  models 
of  composite  electrodes  that  assist  in  predicting  effective  properties 
that  can  be  used  in  larger  scale  continuum  models.  An  anode- 
supported  tubular  SOFC,  such  as  illustrated  in  Fig.  1,  might  be  on 
the  order  of  100  mm  long  and  10  mm  in  diameter.  A  typical  cell 
uses  a  relatively  thick  Ni-YSZ  composite  anode  (around  1  mm)  with 
a  thin  YSZ  dense  electrolyte  (around  10  /xm)  and  a  thin  LSM-YSZ 
composite  cathode  (around  50  /xm).  It  is  well  established  that  SOFC 
performance  depends  greatly  upon  the  structure  of  the  composite 
electrodes  [1].  Thus,  there  is  a  need  to  develop  predictive  models 
that  can  assist  optimizing  the  cell  architecture.  The  present  paper 
is  concerned  with  models  that  describe  solid-phase  ion  and  elec¬ 
tron  transport  at  the  particle  scale  (submicron)  within  composite 
electrodes.  The  particular  objective  is  to  predict  effective  conduc¬ 
tivity  of  a  particle  matrix  as  a  function  of  particle  size  and  packing 
structure.  Although  the  work  is  motivated  primarily  by  solid-oxide 
fuel  cells  (SOFC),  the  approach  and  results  are  considerably  more 
general. 
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Porous  composite  electrodes  involve  three  phases.  The  ion¬ 
conducting  phase  (typically  yttria-stabilized  zirconia,  YSZ,  is 
used  in  both  the  anode  and  cathode)  must  be  a  good  oxygen- 
ion  O2-  conductor,  with  very  little  electronic  conduction.  The 
electron-conducting  phase  (typically  nickel,  Ni,  in  the  anode  and 
strontium-doped  lanthanum  manganate,  LSM,  in  the  cathode)  must 
have  good  electronic  conduction,  with  very  little  ionic  conduc¬ 
tion.  The  pore  phase  must  facilitate  gas  transport.  Typically  each  of 
these  three  phases  occupies  about  one  third  of  the  overall  volume. 
Charge-transfer  reactions  proceed  at  the  three-phase  boundaries 
that  are  formed  at  the  intersections  of  the  three  phases.  The  solid- 
phase  particles  and  the  pore  spaces  are  typically  on  the  order  of  a 
fraction  of  a  micron. 

The  composite  electrode  must  accomplish  several  essential 
functions.  The  electron-conducting  phase  must  provide  a  low- 
resistance  path  from  the  position  where  the  electrons  are  produced 
(or  consumed)  to  the  current  collection.  That  is  to  say  the 
electronic-conducting  phase  must  percolate  through  the  thickness 
of  the  composite  electrode.  Especially  within  a  few  tens  of  microns 
around  the  dense  electrolyte  the  ion-conducting  phase  must  per¬ 
colate,  enabling  ion-conduction  to  the  three-phase  boundaries  [2]. 
When  hydrocarbon  fuels  are  used,  the  Ni  phase  in  the  anode  also 
serves  as  a  reforming  catalyst. 

A  structure  that  facilitates  one  function  can  sometimes  serve  to 
frustrate  another  function.  For  example,  small  particles  and  pore 
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Fig.  1.  Illustration  of  a  tubular  solid-oxide  fuel  cell.  The  balloon  shows  geometrical 
and  functional  aspects  of  the  membrane-electrode  assembly  at  the  particle  scale. 


spaces  tend  to  increase  three-phase-boundary  length,  which  facil¬ 
itates  the  essential  charge-transfer  chemistry.  However,  the  small 
pore  spaces  also  serve  to  increase  gas-transport  resistance,  which 
diminishes  performance.  Small  particles  can  also  decrease  effective 
ionic  and  electronic  conductivity,  which  can  negatively  influence 
performance.  Because  there  are  different  functional  needs  in  dif¬ 
ferent  regions  of  the  composite  electrode,  functional  grading  of 
the  electrode  can  improve  performance.  For  example,  in  the  few 
tens  of  microns  near  the  dense  electrolyte  layer  charge  transfer  is 
most  important  [2].  In  this  region  a  functional  layer  with  fine  parti¬ 
cles  and  small  pore  spaces  is  valuable.  However,  far  from  the  dense 
electrolyte  there  is  no  charge-transfer  chemistry  (i.e.,  all  available 
charge  has  been  transferred  nearer  to  the  dense  electrolyte).  In 
this  region  (usually  the  majority  of  an  anode  support  structure) 
performance  is  improved  by  using  large  pores  that  minimize  gas- 
transport  resistance.  Finding  optimal  electrode  architectures  that 
balance  competing  tradeoffs  can  be  challenging  and  there  is  a  need 
to  understand  quantitatively  the  influence  of  electrode  structure 
on  performance. 

At  the  relatively  large  cell  scale  (e.g.,  the  tube  in  Fig.  1)  mod¬ 
els  cannot  resolve  details  at  the  particle  scale  [3].  The  approach 
in  the  present  paper  begins  by  synthesizing  three-dimensional 
microstructures  using  packed,  partially  overlapping,  spherical 
particles.  From  these  random  arrays  of  particles,  geometrical  char¬ 
acteristics  such  as  coordination  number,  percolation  probability, 
and  three-phase-boundary  length  can  be  evaluated.  With  arrays 
of  submicron  particles  in  a  three-dimensional  domain  of  a  few 
microns,  it  is  possible  model  transport  through  the  fully  resolved 
complex  particle  network.  From  the  solutions  of  such  microscale 
models,  it  is  possible  to  extract  the  effective  conductivities  that  are 
needed  to  model  transport  at  the  macroscale.  The  paper  compares 
effective  properties  derived  from  the  fully  resolved  models  with 
those  derived  from  percolation  theory. 


2.  Prior  literature 

A  significant  and  growing  body  of  literature  is  devoted  to 
understanding  complex  electrode  structure  and  performance  at 
the  microscale.  Experimental  investigations  include  measuring 
three-dimensional  microstructure  using  a  range  of  reconstruction 
techniques.  In  some  cases  the  measurements  are  correlated  to 
electrochemical  performance.  Other  investigators  use  analytical 
and  computation  methods  to  synthesize  and  model  microstruc¬ 
tures. 

Some  investigators  measure  the  structure  of  actual  electrodes. 
Beginning  with  Wilson  et  al.  [4],  several  groups  have  devel¬ 
oped  focused-ion-beam-scanning-electron-microscope  (FIB-SEM) 
experiments  to  reconstruct  composite  electrodes  [5-12].  Izzo  et 
al.  have  developed  similar  electrode  reconstructions  of  actual 
electrodes  using  X-ray  tomography  [13,14].  These  investiga¬ 
tions  provide  a  great  deal  of  quantitative  information  about 
the  complexity  of  actual  electrodes.  Stereological  imaging  and 
reconstruction  provide  high-resolution  measurements  of  three- 
phase-boundary  (TPB)  lengths,  phase  volumes,  specific  surface  and 
interface  areas,  phase  connectivity  (percolation),  and  tortuosity. 
Such  measurements  can  be  used  to  validate  the  simulated  electrode 
microstructures. 

Monte-Carlo  stochastic  simulation  techniques  have  been 
used  to  represent  electrode  microstructures.  Electron-  and  ion¬ 
conducting  particles  can  be  either  randomly  packed  [15-18]  or 
randomly  distributed  on  lattice  structures  [19].  Random  particle 
networks  can  be  further  modified  by  numerical-sintering  algo¬ 
rithms  to  represent  the  composite  electrodes  more  accurately. 
Electrode  structures  can  also  be  generated  by  using  stochastic 
reconstruction  schemes,  which  are  based  upon  spatial  statistical 
correlations  that  are  derived  from  experimental  images  [20,21  ]. 

Reconstructed  microstructures,  either  measured  experimen¬ 
tally  or  synthesized  computationally,  can  be  used  in  models  to 
investigate  the  interactions  of  microstructure  with  the  operating 
performance  [11].  Electrode  reconstructions  can  also  be  param¬ 
eterized  for  incorporation  into  macroscopic-scale  models,  thus 
effectively  bridging  scales  from  microstructure  to  cell  performance 
[10,22].  Important  microstructural  parameters  include  TPB  length, 
phase  volume,  phase  surface  and  interface  areas,  phase  connectiv¬ 
ity  and  conductivity,  and  phase  tortuosity. 

Several  investigations  have  been  concerned  with  evaluating 
phase  tortuosity.  Gostovic  et  al.  [5]  and  Smith  et  al.  [10]  calculated 
tortuosity  by  tracking  the  pore-center  locations  and  evaluating 
effective  path  lengths.  Wilson  et  al.  [4]  evaluated  the  tortuosity  by 
converting  the  measured  three-dimensional  phase  structures  into 
finite-element  meshes  and  solving  steady-state  diffusion  equa¬ 
tions.  An  alternative  approach  is  to  apply  a  lattice  Boltzmann 
method  (LBM)  to  solve  a  diffusion  equation  and  derive  a  measure 
of  phase  tortuosity  [23,20,13,24,12].  Iwai  et  al.  [12]  showed  that 
a  random-walk  method  and  the  LBM  approach  produced  similar 
measures  of  tortuosity. 

Grid-based  counting  techniques  have  been  widely  used  to  esti¬ 
mate  TPB  length  [25,20,6,8,7,26].  Iwai  et  al.  [12]  developed  a 
volume-expansion  method  and  a  centroid  method  to  estimate  the 
TPB  length.  Percolation  theory  and  coordination  number  are  also 
used  to  evaluate  the  active  TPB  length  in  synthesized  microstruc¬ 
tures  [15-18]. 

Although  experimental  reconstructions  provide  specific  details 
of  complex  microstructures,  they  offer  relatively  little  predic¬ 
tive  capability.  Simulation  techniques,  although  not  representing 
actual  microstructures  exactly,  provide  a  valuable  tool  for  pre¬ 
dicting  the  performance  of  alternative  electrode  architectures.  By 
using  measurements  to  validate  microstructural-synthesis  tech¬ 
niques,  the  predictive  ability  of  simulation-based  analysis  is 
improved. 
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3.  Percolation  theory 

At  the  macroscale,  ion  and  electron  transport  is  usually  repre¬ 
sented  using  the  charge-conservation  equations  as 


V-(o*V0m)  +  Km  =  0, 


(1) 


where  a ^  (m  =  ed  for  the  electronic-conducting  phase,  and  m  = 
el  for  the  ionic-conducting  phase)  are  the  effective  conductivities, 
&m  are  the  electric  potentials,  and  Rm  are  the  net  charge-transfer 
rates.  The  challenge  is  to  determine  the  effective  conductivity  in 
the  way  that  represents  the  complexities  of  the  porous  compos¬ 
ite  electrodes.  Effective  conductivities  can  be  approximated  based 
upon  coordination-number  theory  and  the  percolation  theory  as 
[27], 

am  =  CTm[(1  ~  (PsWmPmV ,  (2) 

where  cr°  are  the  conductivities  of  the  bulk  material  (i.e.,  material 
from  which  particles  are  formed),  0g  is  the  porosity,  xj/m  are  the  vol¬ 
ume  fractions  based  on  solids  alone  (e.g.,  ^Ni  =  0.5  and  ^Ysz  =  0.5 
means  that  a  composite  Ni-YSZ  anode  is  composed  of  50%  Ni  and 
50%  YSZ),  and  y  is  called  the  Bruggeman  factor  that  accounts  for 
tortuous  conduction  pathways.  The  Bruggeman  factor  is  typically 
assumed  to  be  y  =  1.5.  The  probability  of  percolation  Pm  repre¬ 
sents  the  probability  that  the  m-phase  particles  form  a  connected 
pathway  through  the  electrode,  enabling  continuous  conduction. 

Bouvard  and  Lange  [28]  evaluated  the  probability  of  an  m-phase 
particle  being  in  a  percolated  cluster  of  the  same  phase  as 


where  the  coordination  number  Zmm  is  the  average  number  of  the 
m-phase  particle  neighbors  of  an  m-phase  particle.  Assuming  all  the 
particles  have  the  same  size,  Zmm  can  be  expressed  as  Zmm  =  \//mZ 
where  Z  is  the  overall  average  coordination  number,  which  is  the 
average  number  of  contacts  that  a  particular  particle  makes  with 
neighboring  particles,  independent  of  the  phase  or  material  identity 
of  the  particles.  It  is  often  assumed  thatZ  =  6  in  a  random  packing 
of  binary  spheres. 

Kuo  and  Gupta  [29]  showed  that  the  percolation  threshold  is 
Zmm  =  1-764,  below  which  no  cluster  connecting  both  ends  of  the 
electrode  is  possible.  Suzuki  and  Oshima  [30]  proposed  a  more 
accurate  expression  for  the  percolation  probability  as, 


Pm 


/ 4.236  -Zmm\2S 
V  2.742  J 


(4) 


which  has  zero  probability  at  the  percolation  threshold. 

The  percolated  three-phase-boundary  (TPB)  length  per  unit  vol¬ 
ume  can  also  be  predicted  from  the  percolation  theory  as  [31  ] 

^TPB  =  ^ed-el^ed^ed-el^ed^eb  (6) 


where  £ed_ei  =  27rr  sin  0  is  the  contact  perimeter  between  an 
ion-conducting  particle  and  an  electron-conducting  particle.  As 
illustrated  in  Fig.  2,  the  particle  radius  is  r  and  the  contact  angle 
is  0.  Zed_e j  is  the  average  coordination  number  between  ion-  and 
electron-conducting  particles,  and  Zed_el  =  ^elZ  when  the  particles 
have  the  same  size.  Ped  and  Pel  are  the  percolation  probabilities  of 
electron-  and  ion-conducting  particles,  respectively.  The  number  of 
the  electron-conducting  particles  per  unit  volume  within  the  entire 
composite  electrode  ned  can  be  evaluated  as  ned  =  ntot^ed  where 
ritot  is  the  total  number  of  particles  per  unit  volume, 


_  1  -  0g 

nt0t  (4/3)77T3  ’ 


(6) 


Fig.  2.  Illustration  of  particle  overlap  angle.  For  different  sized  particles,  the  angle 
0  is  defined  from  the  smaller  particle. 

4.  Fully  resolved  model 

The  approach  used  in  the  present  effort  is  quite  different  from 
percolation  theory.  The  composite  structure  is  synthesized  by 
assembling  a  relatively  small  set  of  particles  (hundreds  of  particles) 
into  a  cubical  region  (typically  a  few  microns  on  a  side).  Neighbor¬ 
ing  particles  may  overlap  to  approximate  sintering.  The  particle 
array  is  then  discretized  (here,  a  tetrahedral  mesh  as  illustrated 
in  Fig.  3)  such  that  a  diffusive  transport  equation  that  is  resolved 
on  the  sub-particle  scale  can  be  solved.  In  this  case,  the  intrinsic 
material  properties  are  used  (i.e.,  the  conductivity  cr° ,  not  an  effec¬ 
tive  conductivity  o ^).  This  approach  accounts  directly  for  particle 
shapes  and  geometrical  effects  such  as  necking  between  particles. 
Because  there  is  a  certain  randomness  in  assembling  the  particle 
arrays,  many  realizations  are  needed  to  deliver  statistically  valid 
results. 

The  process  to  determine  the  effective  conductivity  begins  by 
solving  a  charge-conservation  equation 

V(a°V0m)  =  0  (7) 

within  the  particle  network.  Because  the  intent  here  is  only  to 
determine  effective  conductivities,  the  charge-transfer  chemistry 
is  neglected.  Solution  requires  boundary  conditions  on  all  faces  of 
the  cube.  For  example,  on  the  top  face  3 f2T  the  electrical  potential 
@m  is  imposed  as  @m  =  1  and  on  the  bottom  face  3-f2B,  <Pm  =  0. 
The  side  walls  are  assumed  to  be  symmetry  boundaries,  specified 


Fig.  3.  A  collection  of  randomly  packed  spheres  with  a  tetrahedral  mesh  shown  on 
the  particle  surfaces. 
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as  n  •  V@m  =  0,  where  n  is  a  unit  vector  normal  to  the  face.  With 
a  solution  in  hand,  the  axial  component  of  the  current  through 
the  particle  network  can  be  evaluated  at  either  the  top  or  bottom 
surface  as 

Jfn=  I  ffjn  •  V0mdA  =  [  a^n-V0mdA,  (8) 

J  a^T  J  a^2g 

The  integral  considers  only  the  particle-interface  areas  inter¬ 
secting  the  top  or  bottom  surfaces  through  which  there  can  be 
conduction. 

The  effective  conductivity  is  defined  to  be  the  ratio  of  the  axial 
current  through  the  particle  network  and  the  axial  current  through 
a  solid  cube  of  the  parent  material.  Axial  current  through  a  solid 
cube  is  evaluated  simply  as 


f0  _  /\n o 
Jm  -  t\Om 


d&m 

dx  ’ 


(9) 


where  A  is  cross-sectional  area  of  the  dense  cube.  Thus,  the  effective 
conductivities  are  evaluated  as 


/e 

Jm 

J 0  ' 
Jm 


(10) 


5.  Particle-packing  algorithms 


Fig.  4.  Scanning  electron  microscope  (SEM)  image  of  an  SOFC  membrane-electrode 
assembly  (top)  and  composite  electrode  representations  based  upon  randomly 
packed  spheres. 


Constructing  particle-scale  models  depends  upon  an  algorithm 
to  synthesize  a  particle  and  pore  structure.  For  composite  elec¬ 
trodes,  such  a  structure  is  composed  of  two  particle  types  (i.e., 
ion-conducting  and  electron-conducting)  and  pore  space.  This  is 
most  often  accomplished  using  various  packed-sphere  models 
[32,25],  although  at  least  one  group  uses  rectangular  elements  [11]. 
Monodispersed  spherical  particles  are  used  here,  but  the  approach 
is  fully  capable  of  using  different  shapes  and  polydispersed  sizes. 
In  this  paper  two  different  approaches  are  developed  for  establish¬ 
ing  the  particle  positions;  one  is  based  upon  a  discrete-element 
method  and  the  other  is  a  random  packing. 

Although  random  packing  of  simple  shapes  is  convenient,  it 
should  be  recognized  that  real  structures  may  not  be  so  simply 
represented.  Actual  electrode  fabrication  begins  by  ball  milling 
powders  (nominally  spherical  particles)  together  with  binders.  In 
the  case  of  a  Ni-YSZ  anode,  the  Ni  is  in  the  form  of  nickel  oxide, 
NiO.  The  resulting  slurry  is  formed  into  a  nominal  shape  (e.g.,  tube) 
by  casting  or  extrusion.  The  part  is  then  dried  and  fired  at  high 
temperature,  which  causes  the  particles  to  sinter  together.  In  some 
cases  the  slurry  may  contain  a  pore  former  (e.g.,  graphite),  which 
is  burned  out  during  firing.  Pore  formers  are  used  to  create  pore 
spaces  that  are  larger  than  would  otherwise  be  present  using  the 
functional  powders  alone.  In  a  Ni-YSZ  structure,  the  NiO  must  be 
reduced  to  Ni  before  application  in  the  fuel  cell.  The  NiO  reduction  is 
accomplished  by  reaction  with  H2  at  high  temperature.  Because  the 
NiO  occupies  a  much  larger  volume  than  does  the  Ni,  the  reduction 
process  also  opens  pore  space.  Given  such  a  fabrication  process,  it  is 
easy  anticipate  that  the  resulting  structure  may  not  be  represented 
simply  as  an  assembly  of  packed  spheres. 

Fig.  4  shows  the  comparison  of  a  microscopic  image  of  an  SOFC 
with  a  model  electrode  structures  based  upon  randomly  packed 
spheres.  It  is  evident  that  packed  spheres,  or  any  other  simple 
shape,  are  not  entirely  adequate  in  representing  the  real  struc¬ 
ture.  Nevertheless,  models  based  upon  packed  spheres  do  provide 
valuable  quantitative  insight  about  how  electrode  performance 
depends  upon  microstructure.  In  addition  to  seeking  physically 
representative  packing,  a  further  objective  here  is  to  pack  parti¬ 
cles  in  ways  that  are  consistent  with  assumptions  made  in  models 
based  upon  percolation  theory.  This  provides  a  basis  for  comparing 
resolved-particle  results  with  percolation-theory  predictions. 


5.1.  Discrete  element  method 

The  discrete  element  method  (DEM)  is  a  well-established 
method  for  predicting  the  evolution  of  particle  trajectories, 
accounting  for  interparticle  collisions  and  other  relevant  forces 
[33].  The  method  provides  a  general  framework  for  solving  a  set  of 
ordinary  differential  equations  that  express  the  Newtonian  equa¬ 
tions  of  motion  for  individual  particles.  Martin  et  al.  [34]  first  used 
a  DEM  approach  to  represent  SOFC  anodes,  focusing  on  particle 
deformation  during  sintering. 

Briefly,  as  implemented  here,  the  DEM  packing  algorithm 
involves  the  following  steps: 

•  Define  a  tall  vertical  rectangular  domain,  typically  a  few  microns 
on  short  sides.  That  is  the  height  z  is  around  three  times  greater 
than  the  x  and  y  dimensions. 

•  At  the  top  of  the  rectangular  domain  (x-y  plane)  introduce  small 
(typically  fraction  of  a  micron)  spherical  particles  with  random 
initial  positions  and  velocities.  Each  particle  is  identified  as  either 
an  ion-  or  electron-conducting  particle.  The  particles  fall  under 
the  influences  of  gravity  and  interactions  with  other  particles. 

•  Solve  Lagrangian  equations  of  motion  to  predict  the  particle  tra¬ 
jectory  until  the  particle  comes  to  rest  in  the  lower  portion  of  the 
rectangular  domain. 

•  As  a  particle  moves  its  motion  is  affected  by  the  net  forces  and 
torques  it  encounters  by  colliding  with  other  particles  and  the 
domain  walls. 

•  The  filling  process  is  terminated  when  the  lower  portion  of  the 
rectangular  domain  is  filled  with  particles  at  rest,  forming  a  nom¬ 
inal  cubical  domain. 

•  The  DEM  process  alone  produces  a  set  of  particles  that  touch, 
but  do  not  overlap.  To  simulate  sintering,  the  particle  radii  are 
expanded  to  achieve  a  certain  level  of  overlap. 

For  the  narrow  purposes  of  producing  particle  packing,  the  par¬ 
ticle  contact  forces  are  modeled  using  a  relatively  simple  spring 
and  dashpot  analogy  occurring  over  a  finite  time  interval.  A  more 
detailed  DEM  simulation  could  incorporate  a  more  physically  real¬ 
istic  modeling  of  all  contact  interactions.  However,  because  the 
objective  of  the  present  work  is  to  use  DEM  simply  as  a  means 
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Fig.  5.  Array  of  spheres  produced  from  the  DEM  packing  algorithm.  The  lower  image 
highlights  the  contact  areas  between  the  spheres  after  artificial  sintering. 

to  generate  a  packing  configuration,  there  is  no  need  to  be  con¬ 
cerned  with  potentially  complex  particle-interaction  forces.  The 
particle  masses  as  well  as  the  spring  and  dashpot  parameters  are 
set  somewhat  arbitrarily  to  achieve  reasonable  computation  times. 

The  results  reported  in  this  paper  use  uniform  particles  with 
a  radius  of  0.3  pum.  The  bounding  volume  is  chosen  as  3.3  pum  in 
the  x  and  y  directions.  The  filling  process  continues  until  the  par¬ 
ticles  have  filled  the  domain  up  to  a  z  dimension  about  3.3  |jim. 
Once  the  particle  positions  are  established  the  particle  radii  are 
increased  to  achieve  a  typical  particle-intersection  angle  of  15° 
(Fig.  2).  This  procedure  approximates  a  sintering  process,  and  is 
consistent  with  approximations  that  are  often  made  in  percolation 
theory  [31].  Fig.  5  shows  a  typical  particle-packing  configuration 
produced  by  the  DEM  method.  The  upper  panel  shows  the  parti¬ 
cles  packed  inside  the  cubical  domain.  The  lower  panel  highlights 
the  contact  area  between  particles. 

Porosity,  which  is  an  important  characteristic  of  the  microstruc¬ 
ture,  is  found  by  a  sampling  algorithm.  A  large  sequence  of 
randomly  generated  points  (i.e.,  x,y,  and  z  coordinates)  within  the 
cubical  domain  are  generated.  Each  point  is  evaluated  to  deter¬ 
mine  if  it  is  within  a  solid  particle  or  not.  With  sufficient  numbers 
of  points,  the  porosity  is  established  as  the  ratio  of  points  that  fall 
outside  solid  particles  and  the  points  that  fall  within  solid  particles. 
Points  that  fall  within  one  particle  diameter  of  the  boundaries  are 
discounted  so  as  to  reduce  boundary  effects.  The  average  porosity 
is  found  to  be  about  32%  for  the  DEM  packing  used  here. 

5.2.  Random  particle  packing 

The  second  approach  for  synthesizing  a  microstructure  is  based 
upon  the  random  generation  of  particles  within  a  cubical  domain. 


Valid  configuration 


Invalid  configuration 

Fig.  6.  Valid  and  invalid  particle-packing  configurations. 


Similar  approaches  have  been  used  previously  for  studying  perco¬ 
lating  clusters  [28].  The  algorithm  begins  by  generating  at  random 
a  trial  position  (i.e.,  x,  y,  and  z  coordinates)  for  a  spherical  parti¬ 
cle  within  a  cubical  domain.  The  particle  position  is  tested  against 
a  set  of  constraints  to  determine  its  validity.  If  valid,  the  parti¬ 
cle  is  retained.  Otherwise  the  trial  particle  is  discarded.  In  either 
case,  a  new  particle  position  is  generated  and  tested  for  validity. 
The  process  continues  until  the  cubical  domain  is  filled  with  valid 
particles. 

The  packing  constraints  are  stated  as  follows: 

1 .  A  valid  particle  must  have  contact  or  overlap  with  at  least  three 
neighboring  particles  and/or  the  domain  boundaries. 

2.  The  intersection  perimeter  between  particles  must  be  circular.  In 
other  words,  as  discussed  below,  a  given  particle  cannot  overlap 
simultaneously  with  more  than  one  particle. 

3.  The  intersection  angles  must  be  0°  <  9  <  30°  (Fig.  2),  with  an 
average  overlap  angle  of  15°. 

Fig.  6  shows  examples  of  valid  and  invalid  particle  intersections. 
In  the  upper  panel  sphere  1  satisfies  constraints  with  respect  to 
its  contact  with  spheres  2  and  3  and  surrounding  domain  bound¬ 
aries.  However,  the  configuration  in  the  lower  panel  violates  the 
second  constraint  because  the  intersection  of  spheres  1  and  2  is 
also  intersected  by  sphere  3.  The  second  constraint  is  imposed  to 
generate  configurations  that  are  similar  to  the  ones  obtained  from 
the  DEM-based  method  and  to  be  roughly  consistent  with  com¬ 
mon  assumptions  in  percolation  theory.  However,  this  constraint 
is  not  necessary  and  could  be  relaxed  to  produce  closer  packing  and 
smaller  pore  spaces. 

Filling  a  cubical  domain  using  the  random-packing  algorithm  is 
a  very  lengthy  process  because  most  randomly  generated  positions 
fail  to  satisfy  the  packing  criteria.  However,  the  process  is  fully  auto¬ 
mated.  To  satisfy  the  first  constraint,  stable  particles  are  initially 
positioned  next  to  the  bounding  walls  where  three  or  more  contact 
points  are  easily  found  within  a  nearly  empty  domain.  As  more  and 
more  particles  are  generated,  the  interior  of  the  domain  is  succes¬ 
sively  filled.  There  is  reason  to  expect  from  percolation  theory  that 
particle  overlap  angles  (illustrated  in  Fig.  2)  should  be  around  15° 
[27,35].  However,  because  the  particles  are  generated  randomly,  a 
range  of  interaction  angles  between  0°  and  30°  is  accepted,  yielding 
an  average  overlap  of  about  15°. 

The  porosity  of  the  resulting  microstructure  is  calculated  in 
exactly  the  same  way  as  DEM.  However  due  to  the  various 
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constraints  the  porosities  are  found  to  be  about  47%,  which  is  sig¬ 
nificantly  higher  than  the  porosities  found  via  the  DEM  packing  of 
around  32%.  It  should  be  recognized  that  it  is  possible  to  produce 
lower  porosities  (similar  to  those  with  the  DEM-based  method)  by 
relaxing  the  above  constraints.  However,  doing  so  would  cause  the 
particle-overlap  characteristics  in  the  resulting  microstructure  to 
be  qualitatively  different  from  the  DEM-generated  microstructure. 
Because  it  is  interesting  to  compare  the  two  packing  methods,  the 
random  packing  is  constrained  for  the  purposes  of  this  paper.  A 
porosity  of  47%  is  relatively  high  for  a  typical  SOFC  electrode.  How¬ 
ever,  the  measurements  reported  by  Iwai  et  al.  [  1 2  ]  show  porosities 
in  this  range. 

5.3.  Discretization 


Number  of  realizations 

Fig.  7.  Evolution  of  mean  coordination  number  Z  for  150  and  300  particles  using 
the  DEM  packing. 


Once  the  particle  positions  are  established,  the  next  step  is 
to  discretize  the  structure.  The  particle  coordinates  and  dimen¬ 
sions  of  the  bounding  domain  are  imported  into  the  ANSYS 
Gambit  mesh-generation  software.  Next,  a  set  of  spherical  geomet¬ 
ric  solids  representing  electron-conducting  and  ion-conducting 
particles  are  generated  along  with  the  surrounding  rectangular  vol¬ 
ume.  A  sequence  of  Boolean-intersection  operations  are  performed 
on  these  primitives  producing  the  final  geometry  comprising  of 
overlapping  particles  and  boundaries.  An  unstructured  tetrahedral 
mesh  is  generated  in  the  solid  volumes.  The  complement  of  the 
solid  volumes  formed  by  the  spheres  represent  the  pore  volumes. 
Although  not  used  in  the  present  paper,  the  pore  regions  can  also  be 
meshed  for  solving  the  fluid-transport  equations.  The  three-phase 
boundaries  can  be  automatically  identified  and  the  TPB  lengths 
computed. 

The  unstructured  tetrahedral  meshing  is  entirely  general  and 
is  in  no  way  restricted  to  spherical  particles  or  any  other  array  of 
simple  shapes.  If  a  detailed  geometric  specification  is  available  (e.g., 
from  an  experimental  reconstruction),  the  electrode  structure  can 
be  discretized.  The  simulation  and  analysis  procedure  is  applied  in 
just  the  same  way  as  for  the  spherical-particle  representations. 

6.  Results  and  discussion 

The  primary  objective  of  this  paper  is  to  estimate  effective  con¬ 
ductivity.  However,  the  model  can  also  be  used  to  evaluate  the 
mean  coordination  number  that  is  used  in  percolation  theory  (Eq. 
(4)).  In  percolation  theory  it  is  common  to  assume  an  average 
overall  coordination  number  around  Z  =  6.  Both  the  coordination 
number  and  the  percolation  probability  can  be  explicitly  evaluated 
in  the  present  fully  resolved-particle  model.  Thus  these  two  quanti¬ 
ties  can  be  used  to  provide  further  insight  about  percolation  models 
and  as  partial  validation  of  the  resolved-particle  model. 

6.1.  Coordination  number 

The  coordination  number  is  calculated  directly  once  the  particle 
positions  are  established.  Because  each  instance  of  the  particle- 
generation  process  produces  a  different  geometric  configuration, 
many  such  realizations  are  needed  to  obtain  results  that  are  sta¬ 
tistically  invariant  with  respect  to  the  geometry  and  the  particle 
sample  size.  It  is  expected  that  the  results  obtained  from  simulat¬ 
ing  a  relatively  small  sample  size  for  a  large  number  of  realizations 
should  be  virtually  identical  to  that  obtained  from  simulating  a 
large  sample  size  with  fewer  realizations.  Two  different  sample 
sizes  (150  and  300  particles)  are  simulated  using  the  DEM-based 
model  to  investigate  the  influence  of  particle  sample  sizes  on  the 
coordination  number.  Neither  of  the  packing  algorithms  described 
above  is  limited  by  the  number  of  particles  simulated.  However, 
larger  samples  require  greater  computational  resources  for  mesh¬ 
ing  and  simulation,  limiting  the  number  of  realizations  that  can  be 
practically  considered. 


The  150-particle  samples  are  run  for  600  realizations  and  the 
300-particle  samples  are  run  for  300  different  realizations  result¬ 
ing  in  the  consideration  of  a  total  90,000  particles  for  each  sample 
size.  A  range  of  bins,  initially  empty  and  representing  discrete 
coordination  numbers,  are  set  up  at  the  beginning  of  the  computa¬ 
tions.  Following  a  given  instance  of  the  particle-generation  process, 
the  coordination  number  for  each  particle  is  computed  and  the 
appropriate  bin  count  is  incremented.  Next,  the  mean  coordina¬ 
tion  number  is  computed  for  all  particles  and  a  running  average 
is  maintained  over  all  realizations.  The  process  is  repeated  until  a 
quasi-steady  value  of  the  mean  coordination  number  is  achieved. 
This  procedure  effectively  determines  the  number  of  realizations 
(and  hence  particle  sample  size)  required  to  achieve  a  statistically 
invariant  value  of  the  mean  coordination  number.  Additionally  the 
probability  distribution  function  (PDF)  of  the  coordination  number 
is  obtained  as  a  result  of  the  simulations.  Only  particles  which  are 
situated  at  least  one  particle  diameter  away  from  the  walls  are  con¬ 
sidered  in  the  coordination-number  calculation,  thus  minimizing 
wall  effects. 

Fig.  7  shows  the  typical  evolution  of  the  mean  coordination 
number  as  a  function  of  the  number  of  realizations  for  150  and 
300  particles  with  the  DEM  packing.  The  mean  coordination  num¬ 
ber  asymptotically  approaches  Z  ^  6.7.  This  result  confirms  that 
the  strategy  of  performing  a  large  number  of  simulations  with  a 
limited  sample  size  is  theoretically  sound.  The  final  value  of  the 
coordination  number  is  consistent  with  results  reported  by  others 
[28]. 

Fig.  8  shows  comparative  histograms  produced  by  the  150-  and 
300-particle  simulations  with  the  DEM  packing.  Once  again  it  is  evi¬ 
dent  that  the  PDF  of  coordination  numbers  obtained  from  the  two 
sample  sizes  are  virtually  identical.  Similar  trends  are  also  observed 
for  the  random-packing  algorithm.  Based  upon  experience  with  the 
DEM-packing  and  the  random-packing  algorithms,  all  subsequent 
simulations  were  performed  with  a  sample  size  of  1 50  particles  and 
600  realizations. 
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Fig.  8.  Probability  of  mean  coordination  number  Z,  comparing  samples  of  150  and 
300  particles  with  DEM  packing. 
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Fig.  9.  A  typical  percolating  cluster. 


6.2.  Percolation  probability 

Percolation  probability  plays  an  important  role  in  evaluating 
effective  three-phase-boundary  length  and  effective  conductiv¬ 
ity.  Once  a  set  of  particles  have  been  generated,  clusters  can  be 
identified,  and  the  percolation  probability  evaluated.  Comparing 
the  directly  evaluated  percolation  probabilities  with  percolation 
theory  provides  some  partial  validation  of  the  resolved-particle 
algorithms. 

An  algorithm  is  formulated  to  identify  distinct  clusters  of 
contiguous  particles  of  a  given  type  (i.e.,  ion-conducting  or 
electron-conducting).  List-based  data  structures  are  established  to 
represent  contiguous  particles  or  clusters.  Beginning  with  a  par¬ 
ticular  particle,  the  algorithm  keeps  track  of  all  particles  that  are 
connected  directly  or  indirectly,  thereby  forming  a  contiguous  clus¬ 
ter.  The  process  is  repeated  for  all  particles.  Finally  all  particles  in 
the  domain  are  assigned  to  unique  clusters.  Once  individual  clus¬ 
ters  are  identified,  they  are  compared  with  the  domain  extents  to 
isolate  only  those  clusters  that  percolate  throughout  the  domain 
in  either  the  x-,  y-,  or  z-direction.  Fig.  9  shows  a  typical  cluster 
spanning  the  domain  in  the  x-direction,  thus  forming  a  continuous 
conduction  pathway.  The  percolation  probability  for  a  given  phase 
volume  fraction  is  calculated  as  the  fraction  of  cases  when  one  or 
more  percolating  clusters  are  found. 

Fig.  10  shows  the  percolation  probabilities  as  functions  of  phase 
volume  fraction.  The  resolved-particle  models  predict  that  the  per¬ 
colation  probability  increases  more  slowly  as  a  function  of  phase 
volume  fraction  than  the  percolation  theory  (Eq.  (4)).  Moreover, 
the  particle  models  predict  higher  percolation  probability  at  low 
phase  volume  fraction.  In  other  words,  the  percolation  threshold 
occurs  at  lower  phase  volume  fraction.  Nevertheless,  the  resolved 
models  predict  the  rapid  increase  in  percolation  probability  in  the 
volume-fraction  range  0.3  <  \l/m  <  0.5,  indicating  that  the  com- 


Fig.  10.  Percolation  probability  as  functions  of  solid-phase  volume  fraction,  \[fm. 
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Fig.  11.  Specific  three-phase-boundary  (TPB)  lengths  A,^pB  predicted  by  percolation 
theory  (i.e.,  Eq.  (5))  compared  to  the  resolved-particle  model  with  DEM  packing.  The 
percolated  TPB  lengths  predicted  by  percolation  theory  use  two  alternatives  (i.e.,  Eq. 
(3)  or  (4))  to  evaluate  the  percolation  probabilities. 

putational  approach  provides  a  reasonable  representation  of  the 
microstructure. 

6.3.  Three-phase-boundary  length 

Accurate  prediction  of  the  three-phase-boundary  (TPB)  length 
is  an  important  aspect  of  modeling  charge-transfer  chemistry  in 
SOFCs.  The  TPB  length  is  evaluated  by  discretizing  the  circles  that 
are  formed  around  the  intersections  between  an  ion-conducting 
and  electron-conducting  particles,  and  then  counting  the  number 
of  points  on  that  boundary  that  are  also  shared  by  the  adjoining 
pore  volume. 

Fig.  11  shows  TPB  lengths  (i.e.,  length  per  unit  volume)  pre¬ 
dicted  by  percolation  theory  and  the  DEM  packing  algorithm.  For 
a  three-phase  boundary  to  be  effective  in  charge  transfer,  both 
ions  and  electrons  must  be  available.  That  is,  there  must  be  per¬ 
colation  for  both  the  ion-  and  electron-conducting  phases  to  the 
TPB.  The  upper  panel  in  Fig.  1 1  is  concerned  with  the  entire  TPB 
lengths,  independent  of  whether  there  is  percolation  or  not.  The 
resolved  model  and  the  percolation  theory  (i.e.,  Eq.  (5),  assuming 
Pe i  =Pe d  =  1 )  agree  well,  indicating  that  the  geometric  evaluation 
of  complete  TPB  length  is  consistent  with  percolation  theory.  The 
lower  panel  of  Fig.  1 1  is  concerned  with  the  percolated  TPB  length. 
Two  percolation- theory  results  are  shown,  using  Eqs.  (3)  (Bouvard 
and  Lange  [28])  and  (4)  (Suzuki  and  Oshima  [30])  to  evaluate  the 
percolation  probabilities.  Owing  to  lack  of  percolation  for  one  of  the 
phases  at  high  and  low  volume  fractions,  the  percolated  TPB  shows 
a  much  narrower  distribution  than  the  complete  TPB  length.  The 
DEM  simulations  agree  reasonably  well  with  the  percolation  the¬ 
ory,  although  the  distribution  functions  are  wider  and  the  peak  TPB 
length  is  somewhat  lower. 

Fig.  12  shows  TPB  lengths  predicted  by  the  random-packing 
algorithm.  Because  the  porosity  is  higher  than  in  the  DEM  simula¬ 
tions,  the  TPB  lengths  are  generally  smaller  for  the  random  packing. 
The  widths  of  the  percolated  TPB  distributions  are  comparable  for 
the  DEM  and  random-packing  cases,  but  the  predicted  peak  per¬ 
colated  TPB  lengths  are  considerably  lower  than  the  percolation 
theory  for  the  random-packing  simulations. 
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Fig.  12.  Specific  three-phase-boundary  (TPB)  lengths  A.^pB  predicted  by  percolation 
theory  (i.e.,  Eq.  (5))  compared  to  the  resolved-particle  model  with  random  packing. 
The  percolated  TPB  lengths  predicted  by  percolation  theory  use  two  alternatives 
(i.e.,  Eq.  (3)  or  (4))  to  evaluate  the  percolation  probabilities. 

6.4.  Effective  conductivity 

The  resolved-particle  models  evaluate  the  effective  conduc¬ 
tivity  directly.  However,  based  upon  percolation  theory,  one  can 
anticipate  certain  functional  dependencies.  Eq.  (2)  shows  that 


Fig.  13.  Conduction  in  two  different  percolating  clusters.  Colors  indicate  electric 
potential,  with  red  being  high  and  blue  being  low.  (For  interpretation  of  the  refer¬ 
ences  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the 
article.) 


Fig.  14.  Predicted  effective  conductivities  cr^  as  functions  of  solid-phase  volume 
fraction. 

the  effective  conductivity  depends  upon  the  percolation  prob¬ 
ability  and  the  tortuosity  of  the  conduction  pathways  (i.e.,  the 
Bruggeman  factor,  y).  The  probability  of  a  percolating  cluster  of 
a  certain  particle  type  increases  with  increasing  particle  volume 
fraction. 

The  resolved-particle  models  are  applied  for  a  wide  range  of 
phase  volume  fractions.  Once  a  percolating  cluster  is  identified,  a 
pure  diffusion  problem  is  solved  for  a  generic  scalar  (e.g.,  electric 
potential)  with  Dirichlet  boundary  conditions  specified  at  each  end 
(Eq.  (7)).  Fig.  13  shows  typical  surface  distributions  of  the  scalar 
in  clusters  for  two  different  realizations  (particle  configurations). 
To  obtain  statistically  invariant  results,  the  simulations  are  run  for 
twenty  distinct  realizations  corresponding  to  each  phase  volume 
fraction.  The  effective  conductivity  is  defined  as  the  ratio  of  the  total 
current  through  the  conducting  clusters  for  a  given  phase  volume 
fraction  and  the  ideal  current  through  the  cubical  volume  (Eq.  ( 1 0)). 

Fig.  14  shows  predicted  effective  conductivities  as  functions  of 
phase  volume  fraction  for  both  DEM  and  random  packing.  The  DEM 
packing  produces  average  porosity  of  0g  =  0.32  and  the  random 
packing  produces  average  porosity  of  </>g  =  0.47.  The  figure  also 
shows  comparisons  with  “standard”  percolation  theory  (i.e.,  Eq. 
(2)),  but  using  different  values  of  the  Bruggeman  factor.  It  is  evident 
that  the  commonly  used  value  of  y  =  1.5  produces  much  higher 
effective  conductivities  than  those  predicted  by  the  fully  resolved 
models.  A  value  of  y  «  3.5  is  consistent  with  the  resolved-particle 
simulations. 

There  is  some  evidence  that  the  effective  conductivity  predicted 
by  conventional  percolation  theory  (i.e.,  Bruggeman  exponent  of 
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y  =  1.5)  may  be  too  high.  Recent  publications  by  Iwai  et  al.  [  1 2  ]  and 
Wang  et  al.  [36]  are  consistent  with  the  present  resolved-particle 
findings,  showing  significant  deviations  from  percolation  theory. 
Iwai  et  al.  propose  that  the  effective  conductivity  may  be  evaluated 
as  [12] 


_0  (1  —  0g  Wm 


(11) 


where  rm  is  a  tortuosity  factor.  In  this  equation  (1  -  0g)/^m  is  the 
fraction  of  the  entire  volume  occupied  by  phase  m,  whereas  \l/m  is 
the  volume  fraction  of  the  solid  phases  alone.  Although  the  tortu¬ 
osity  factors  reported  by  Iwai  et  al.  are  strongly  anisotropic,  even 
the  least  tortuous  paths  produce  effective  conductivities  that  are 
much  lower  than  those  predicted  by  Eq.  (2). 

Wang  et  al.  [36]  formulated  stochastic  representations  of  PEMFC 
composite  cathodes,  assuming  the  microstructure  to  be  repre¬ 
sented  by  randomly  generated  cubes.  They  simulated  transport 
and  charge  transfer  using  a  fully  resolved  model  similar  to  the 
one  developed  in  the  present  paper.  Comparing  the  fully  resolved 
results  to  a  one-dimensional  macromodel  suggested  a  Bruggeman 
factor  between  3.5  and  4.5.  Thus,  the  Wang  et  al.  results  are  also 
consistent  with  the  present  models. 

The  Bruggeman  factor  y  is  intended  to  account  for  two  separate 
phenomena.  One  is  the  convoluted  nature  of  conduction  pathways 
and  second  is  the  possibility  of  isolated  non-percolating  clusters 
that  do  not  contribute  to  conduction.  Although  the  percolation 
probability  (i.e.,  the  probability  that  at  least  one  cluster  will  perco¬ 
late)  is  unity  for  >  0.5,  there  may  be  islands  of  non-percolating 
particles  that  are  not  taken  into  account  by  Eq.  (2 ).  The  fully  resolved 
models  capture  such  islands  and  thus  predict  significantly  lower 
values  of  effective  conductivity  than  the  percolation  theory. 


7.  Conclusions 


A  fully  resolved-particle-based  model  of  composite  electrodes 
has  been  developed  and  applied  to  predict  average  coordination 
numbers,  percolation  probabilities,  TPB  lengths,  and  effective  con¬ 
ductivities  of  the  packed-particle  networks.  Two  different  packing 
algorithms  are  used  to  synthesize  microstructures  from  spheri¬ 
cal  particles,  with  both  algorithms  producing  similar  results.  Once 
a  particle-based  microstructure  is  generated,  the  entire  particle 
matrix  is  discretized  with  a  tetrahedral  mesh  network.  A  charge- 
conservation  equation  is  solved  to  predict  current  through  the 
particle  network.  These  results  are  then  used  derive  an  effective 
conductivity  that  is  represented  as  a  fraction  of  the  intrinsic  mate¬ 
rial  conductivity  and  microstructural  characteristics  of  the  particle 
network. 

The  particle-packing  algorithms  are  random,  causing  each 
realization  of  a  particle  packing  to  be  different.  Thus,  multiple  sim¬ 
ulations  are  needed  to  develop  statistically  valid  results.  The  results 
show  that  a  few  hundred  particles  with  a  few  hundred  realizations 
are  sufficient  to  reduce  statistical  variations  to  within  acceptable 
levels. 

Coordination  numbers,  percolation  probabilities,  and  TPB 
lengths  can  be  evaluated  based  upon  packing  geometry  alone. 
These  results  are  found  to  be  consistent  with  percolation  theory, 
serving  as  partial  validation  of  the  model.  The  effective  conduc¬ 
tivities  that  are  derived  from  the  fully  resolved  simulations  are 
especially  valuable  in  larger  scale  models  where  the  particle  scale 
cannot  be  resolved.  The  effective  conductivities  are  found  to  be 
significantly  smaller  than  those  predicted  by  standard  percolation 
theory.  However,  using  a  Bruggeman  factor  of  y  =  3.5  instead  of 
the  commonly  used  y  =  1.5  brings  the  percolation  prediction  in 
line  with  the  fully  resolved  results. 


The  results  in  the  present  paper  are  concerned  with  uniform- 
diameter  spherical  particles.  One  reason  for  this  choice  is  to 
compare  results  with  percolation  theory.  However,  the  resolved- 
particle  modeling  approach  accommodates  different  particle 
shapes  and  sizes,  thus  providing  the  capability  to  explore  and  eval¬ 
uate  innovative  electrode  architectures.  Moreover,  the  model  can 
be  extended  easily  to  evaluate  fluid-phase  transport  through  the 
pore  volumes. 
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